%clear all;

% load and plot data
load hr1.mat;
data_length = length(rr);
[n,d]=size(rr); %numero punti e dimensione del punto, qui 1D

%determino i range per il plotting
minr=1;
maxr=max(rr);
range =linspace(minr,maxr,500)';

figure(1);

% plot hrv
subplot(1,5,1);
plot(rr, 'b');
title('HRV R-R data')

subplot(1,5,2);
bins=100; % questi sono i bin
plot(hist(rr,bins));
title('R-R histogram');



% define mixture model
ncentres = 3; input_dim = 1;
mix = gmm(input_dim,ncentres,'spherical');
options = foptions;
%options(4) = 10;
options(14) = 5;	% Just use 5 iterations of k-means in initialisation

% initialise mixture from data
%mix = gmminit(mix,rr,options);

% em happiness
maxIter = 30; options(3) = 0.1; options(14) = maxIter;
%[mix2, options, errlog] =  gmmem(mix, rr, options);


% plotto la mistura approssimata
subplot(1,5,3);
title('Mistura di Gaussiane stimata');
hold on


 ymix=zeros(size(range));
 for k=1:3
     y=zeros(size(range));
     
     g = normpdf(range, mix2.centres(k), mix2.covars(k));
     
     %g = gauss(mix2.centres(k), mix2.covars(k), range)
     y= mix.priors(k)*g;
     y=g;
     ymix=ymix+y
     plot(range,y,'b-'), grid on;
     hold on,
 end
plot(range,ymix,'r-'), grid on;
hold off

%calcolo la probabilit�  dei dati secondo il modello a mistura
% P(rr)=\sum_k P(k) P(rr; mu_k, sigma_k)
P_rr=gmmprob(mix2, range);


subplot(1,5,4);
plot(P_rr)

subplot(1,5,5);
title(' GMM stimata vs distribuzione empirica');
hold on
%normalizzo tutto tra zero e uno per plottare e confrontare con l'empirica
maxP_rr=max(P_rr); minP_rr=min(P_rr);
P_rrplot=(P_rr-minP_rr)./(maxP_rr-minP_rr);


[a,b]=hist(rr,bins); 
freq = a./sum(a); %CALCOLO LA DISTRIBUZIONE EMPIRICA in [0,1] per poterla confrontare con Y

maxF=max(freq); minF=min(freq);
Fplot=(freq-minF)./(maxF-minF);


plot(b,Fplot,'ro'); %plotto empirica
plot(range,P_rrplot,'k-'); %plotto stimata




% % plot results
% figure(2);
% subplot(1,2,1);
% g1 = normpdf(1:max(rr), mix2.centres(1), mix2.covars(1));
% g2 = normpdf(1:max(rr), mix2.centres(2), mix2.covars(2));
% g3 = normpdf(1:max(rr), mix2.centres(3), mix2.covars(3));
% hold on
% 
% plot(g1, 'b');
% plot(g2, 'r');
% plot(g3, 'g');
% title('plot gaussiane con dati lineari')
% hold off
% 
% % ????????????
% subplot(1,2,2);
% g1 = normpdf(rr, mix2.centres(1), mix2.covars(1));
% g2 = normpdf(rr, mix2.centres(2), mix2.covars(2));
% g3 = normpdf(rr, mix2.centres(3), mix2.covars(3));
% hold on
% plot(g1, 'b');
% plot(g2, 'r');
% plot(g3, 'g');
% title('plot gaussiane con dati rr')
% 
% hold off

